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We use molecular dynamics simulations to investigate translational and rotational diffusion in 
a rigid three-site model of the fragile glass former ort/io-terphenyl, at 260 K < T < 346 K and 
ambient pressure. An Einstein formulation of rotational motion is presented, which supplements 
the commonly-used Debye model. The latter is shown to break down at supercooled temperatures 
as the mechanism of molecular reorientation changes from small random steps to large infrequent 
orientational jumps. We find that the model system exhibits non-Gaussian behavior in translational 
and rotational motion, which strengthens upon supercooling. Examination of particle mobility re- 
veals spatially heterogeneous dynamics in translation and rotation, with a strong spatial correlation 
between translationally and rotationally mobile particles. Application of the Einstein formalism to 
the analysis of translation-rotation decoupling results in a trend opposite to that seen in conven- 
tional approaches based on the Debye formalism, namely an enhancement in the effective rate of 
rotational motion relative to translation upon supercooling. 



I. INTRODUCTION 

Molecular motion in supercooled liquids has received 
much scrutiny from experiments and simulations in re- 
cent years 0, . These studies have uncovered a rich 
phenomenology not present at temperatures above the 
melting point. One distinguishing feature is commonly 
known as dynamic heterogeneity. This refers to the 
presence of transient spatially separated regions with 
vastly different relaxation times EL Observed in experi- 
ments HHH and simulations [1 SHU 113, these do- 
mains are separated by only a few nanometers, but can 
differ by up to five orders of magnitude in their rate of 
relaxation Closely related to dynamic heterogeneity 
is a non-Gaussian distribution of particle displacements 
at times intermediate between the ballistic and diffusive 
regimes of molecular motion 8] . This behavior is a com- 
mon aspect of supercooled liquids and has been used ex- 
tensively to detect dynamic heterogeneity in computer 
simulations [|| Q • Dynamic heterogeneity has also been 
invoked to explain the decoupling between translational 
diffusion and viscosity, and between rotational and trans- 
lational diffusion in deeply supercooled liquids [ll ITll Il2| | . 
At T > l-2T g , where T g is the glass transition temper- 
ature, the translational diffusion coefficient, D t , and the 
rotational diffusion coefficient, D r , are proportional to 
Tr\~ x , where r\ is the shear viscosity. This is in accord 
with the Stokes-Einstein (SE) equation for translational 
diffusion 
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and its rotational counterpart, the Debye-Stokes-Einstein 
(DSE) relation 
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In both expressions k B is Boltzmann's constant and R is 
an effective hydrodynamic radius of the diffusing parti- 
cle. The fact that these equations, which describe diffu- 
sion of a Brownian particle in a hydrodynamic continuum 
with viscosity rj, apply at the molecular scale is remark- 
able. Yet, as liquids enter the deeply supercooled regime, 
T < 1.2T g , D t exhibits a weaker temperature dependence 
while D r continues to adhere to the DSE relation. Exper- 
iments on ori/io-terphenyl (1,2-diphenylbenzene, OTP) 
reveal that D t is approximately proportional to Tr)~ a s in 
this regime and that the SE equation underpredicts D* by 
as much as two orders of magnitude at T g + 3 K [Hi. Il4j . 

The goal of the present work is to investigate numer- 
ically the diffusive phenomena in a model supercooled 
liquid. Special emphasis is placed on rotational motion, 
an important aspect of supercooled liquid dynamics that 
has received comparatively less attention than its trans- 
lational counterpart in existing investigations of dynamic 
heterogeneity and non-Gaussian behavior. Examples of 
recent studies of rotational dynamics in supercooled liq- 
uids include refs. 0,0,0,0,0- Inclusion of rota- 
tional degrees of freedom provides the means to study 
computationally translation-rotation decoupling, a phe- 
nomenon observed experimentally in many fragile glass 
formers. To investigate these topics we perform molecu- 
lar dynamics simulations of the rigid three-site Lewis and 
Wahnstrom model of OTP [24( at temperatures spanning 
the warm thermodynamically stable liquid to deeply su- 
percooled states. 

This paper is organized as follows. Section ILTl provides 
the background on the two formalisms used here to char- 
acterize rotational diffusion; Section IlIII extends metrics 
for dynamic heterogeneity and non-Gaussian behavior of 
translational motion to rotation. Section llVl presents the 
results obtained from molecular dynamics simulations of 
OTP and discusses their significance. The major con- 
clusions arrived at in this work and the open questions 
arising as a result of our study are listed in Section Ivl 
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II. ROTATIONAL FORMALISM 

The most common framework for exploring rotational 
motion originates with the work of Debye p0| . The un- 
derlying physical picture views rotational diffusion as a 
succession of small, random processes. A unit vector 
fixed to the center of mass of a rotating molecule would 
then undergo a random walk on the surface of a sphere. 
Solving the differential equation for the evolution of the 
probability P(ip, t) that a molecule experiences a net an- 
gular displacement tp during a time t then yields |2l| 

m*) = E (^WosV^le-^ 1 )^ (3) 
1=1 ^ ' 

Here P; is the Z th Legendre polynomial and D r is the ro- 
tational diffusion coefficient, with units of inverse time. 
The angular displacement is defined as ip(t) — cos -1 [u(t) ■ 
u(0)], where u is the unit vector fixed in the molecu- 
lar frame. The first two Legendre polynomials, I = 1 
and 1 = 2, can be related to several experimental tech- 
niques including infrared absorption, Raman scattering, 
and NMR 22| . These experiments often report rotational 
correlation times, t; = [1(1 + l)D r ] _1 , in place of a dif- 
fusion coefficient. These are straightforwardly calculated 
from 

/>OC 
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At low temperatures, Pi[cosijj(t)] decays slowly and long 
simulations are needed to reach the time necessary to 
accurately determine r; from equation 0] One may de- 
rive an equivalent relation from equation [3] and obtain a 
rotational correlation time from 

rr 1 = ~MPi[<x*m]) (5) 

in the region where \n(Pi[cosip(t)]) is linear in time. We 
have verified that equations E| and give similar values of 
t\ at high temperature where P [cos ?/>(£)] decays quickly. 
We refer to this formulation of rotational motion as the 
Debye model. 

The Debye model is well-suited for examining the 
rotation of dipolcs and diatomic or other rigid linear 
molecules because these systems provide a natural choice 
for the unit vector u. However, in molecules with more 
than two rotational degrees of freedom (any rigid non- 
linear molecule), at least two orthogonal unit vectors 
are required for a full description of rotational motion. 
As will be shown below, rotational motion along differ- 
ent directions can exhibit widely differing characteristics 
at low enough temperatures. Additionally, the Debye 
model has limited utility in examining dynamic hetero- 
geneity and non-Gaussian behavior because the angu- 
lar displacement, ip, is bounded between and 7T. As 
will be shown in Section II I II not only is an unbounded 
displacement critical to studying dynamic heterogeneity, 



but rotational motion in the deeply supercooled regime 
deviates appreciably, at least for the Lewis and Wahn- 
strom model considered here, from the physical picture 
of small uncorrelated angular displacements underlying 
the Debye model. In light of these facts, the following 
alternative rotational formalism is introduced. An angu- 
lar displacement may be obtained in a manner analogous 
to that of translational displacement by integrating the 
angular velocity vector [23| . 

A0i(t) = <pi(t) - 0i(Q) = f t3i(t') dt' (6) 

Jo 

Here, uji and Aipi are the angular velocity and angu- 
lar displacement vectors for particle i respectively. Note 
that Atpi is unbounded. For a generic rigid body with 
three rotational degrees of freedom each component de- 
scribes rotation about a specific principal axis fixed in the 
molecular frame and originating at the center of mass. 
The components are denoted by: uji = [u>f, u>^, oj?] and 

In this approach, rotational motion is anisotropic; each 
direction describes a specific rotational movement of the 
molecule. Associated with each rotational direction is a 
principal moment of inertia. These in general are not 
equivalent, and it is therefore appropriate to define diffu- 
sion coefficients for each degree of freedom, rather than 
for the molecule as a whole. As a starting point, con- 
sider the Einstein relation for the translational diffusion 
coefficient, D t , in one dimension. 

1 N 

A=li^E<[^«] 2 > (7) 

i=l 

Here Axi(t) is the displacement of particle i in the re- 
direction at time t, and the sum is over all molecules. A 
rotational diffusion coefficient in the ^"-direction, P", 
can be defined analogously as 

1 N 

^-^^E([A<W] 2 > (8) 

i=l 

where a £ [x, y, z]. We refer to this description of rota- 
tional motion as the Einstein formulation, but note that 
in the rotational case the directions ip x , tp y , and tp z are 
in the molecular frame. An interesting feature of this 
formulation is that by appropriately selecting the orien- 
tational axes fixed in the molecular frame, the different 
diffusion coefficients will correspond to specific changes 
in orientation of the molecule. The axes chosen in this 
paper are illustrated in Figure ^ 

Each of the above models of rotational motion has 
its advantages. The Debye model forms the basis of 
several experimental techniques (primarily via measure- 
ments of P2[cosip(t)]), and has therefore played an im- 
portant role in studies of translation-rotation decou- 
pling miii mm. In contrast, the Einstein formula- 
tion has only been used in simulations p^UrR I23I I2H |28| . 
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FIG. 1: Illustration of rotational directions in the Lewis and 
Wahnstrom model of OTP. In this model each phenyl ring is 
represented by a Lennard- Jones site, and the three sites con- 
stitute a rigid isosceles triangle with a vertex angle of 75° |24j . 
The unit vector (u) used here to study rotation according to 
the Debye model is shown. 



but as will be shown below, it offers important ad- 
vantages for the numerical investigation of dynamics in 
deeply supercooled liquids, where the physical basis of 
the Debye approximation appears to break down. 



III. NON-GAUSSIAN METRICS AND 
DYNAMIC HETEROGENEITY 

The notion of dynamic heterogeneity implies the exis- 
tence of a time scale intermediate between the ballistic 
and diffusive regimes over which particles of like mobil- 
ity are spatially correlated. Therefore, in order to study 
dynamic heterogeneity one must quantify "mobility" and 
identify this time scale. Most studies have defined mobil- 
ity starting from the individual particle displacement vec- 
tor Afi(At) and its associated scalar quantity Afj(At). 
To facilitate comparison between translational and rota- 
tional diffusion, and to connect with a recently developed 
diffusion formalism j29j, we elect to utilize the Carte- 
sian component of the displacement vector, Axi(At) for 
translation (the y and z directions are equivalent and are 
included as independent samples in all our calculations), 
and the individual angular displacements Aipf(At), for 
rotation, as the quantities of interest. 

The work of Kob et al. revealed that the dynamic 
heterogeneity of translational motion observed in com- 
puter simulations of a binary glass-former is associated 
with a non-Gaussian probability distribution of particle 
displacements. For short times displacement may be ap- 
proximated by Axi « vf At resulting in a probability dis- 
tribution that is Gaussian due to the Maxwell-Boltzmann 
distribution of velocities. 
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Also, at sufficiently long times diffusion is purely random 
and the displacement distribution asymptotically adopts 



a Gaussian form. 
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At intermediate times where these approximations are 
invalid, the behavior of P(Axi) has been shown to be 
substantially non-Gaussian, reaching a maximum devi- 
ation from the Gaussian form at a time to be denoted 
by At*. This time serves as the appropriate time scale 
on which to study dynamic heterogeneity of translational 
motion, and corresponds approximately to the beginning 
of the long-time diffusive regime. 

The extension of this idea to rotational motion is 
straightforward. At short times angular displacement 
may be approximated by A<pf ss wf At and the proba- 
bility distribution of angular displacements is accordingly 
Gaussian due to the Maxwell-Boltzmann distribution of 
angular velocities. 



P(Acpf) 
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Here I a is the moment of inertia for a rotation in the 
(/? Q -direction with a 6 [x, y, z\. In analogy to translation, 
the distribution of angular displacements also converges 
to a Gaussian form at long times. 
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Similarly to translation, we will show that at times inter- 
mediate between these approximations P(Aipf) is non- 
Gaussian, and we will focus on the time of maximum 
non-Gaussian behavior At*. We use this time to exam- 
ine dynamic heterogeneity of rotational motion. The no- 
tation At* is used here generically, and we will see that 
At* for translation is not in general the same as At* for 
rotation. 

Determination of At* requires an appropriate metric 
of non-Gaussian behavior. This is commonly done by 
using the ratio of the second and fourth moments of the 
appropriate probability distribution of particle displace- 
ment 30] . A useful non-Gaussian parameter is 
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for translation, and 



al (At) = 
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for rotation. Each parameter is defined such that Gaus- 
sian behavior yields a\ = ol t 2 = 0. This can be easily 
verified from equations an d El The time at which 
a\ and a? reach their maximum value identifies At* for 
translation and rotation, respectively. 

Having specified a time scale, At*, particles may be 
classified as "mobile" and "immobile" using an appropri- 
ate displacement cutoff. Objective specification of this 
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cutoff is challenging and previous studies have explored 
several options. We elect to use a method developed by 
Shell et al. @ which determines a displacement cutoff 
that can be easily adapted to the rotational formalism 
used here. This technique fits a two-Gaussian function 
of the form 

P(z) = fG(z;* 1 ) + (l-f)G(z;<r 2 ) (15) 

to the probability distribution of displacements at the 
time At*. Here / is a weighting parameter with < 
/ < 1) G(z:a) is a Gaussian distribution in z with zero 
mean and standard deviation a, and z signifies Ax or 
Aip a . This functional form reproduces the measured 
distribution at At* well with only a slight underesti- 
mation of the tails. A critical displacement value can 
then be defined from the standard deviations of the 
two-Gaussian fit. For translation, a reasonable three- 
dimensional cutoff is Ar * = V3(ai + ct 2 )/2 m Th e 
\/3 factor accounts for calculation of the standard devi- 
ations from a one-dimensional probability distribution, 
since ([Ar] 2 ) = 3([A.t] 2 ). For each rotational direc- 
tion, a one-dimensional cutoff is defined as A(p a * — 
(a i + a- 2 )/2. Particles are classified as translationally 
mobile if Ari(At*) > Ar* and rotationally mobile in the 
^-direction if A<pf (At*) > A<p a * . 

IV. SIMULATION RESULTS AND DISCUSSION 

We have performed molecular dynamics simulations 
of the Lewis and Wahnstrom model for OTP [24|. In 
this model each benzene ring is represented by a sin- 
gle Lennard- Jones (LJ) site on a rigid isosceles trian- 
gle. The two short sides of the triangle are one LJ 
diameter, a, in length (0.483 nm), and the long side 
is 1.217 LJ diameters long (0.588 nm), giving a vertex 
angle of 75°. The sites on different molecules inter- 
act pairwise-additively with a LJ interaction energy of 
e/^s = 600 K. The chosen directions of the molecular 
frame give I x : I y : I z = 1 : 1.77 : 2.77. Calculation 
runs are performed in the NVE ensemble for N = 324 
molecules (972 LJ sites) at a range of temperatures of 
260 K < T < 346 K. The density at each temperature 
is given in Table [I] and corresponds to the equilibrium 
densit y fo r the model at 1 bar as determined by Rinaldi 
et al. 31] . With the exception of quantities reported in 
Table U all values are reported in reduced units with the 
mass of an OTP molecule, characteristic interaction en- 
ergy e, and atomic diameter a as the reference units. The 
reduced time unit is then 3.19 ps. The rigid body equa- 
tions of motion were integrated using an iterative quater- 
nion algorithm 32] with a step size of 0.001 reduced time 
units. To improve statistics, multiple time origins sepa- 
rated by 0.1-0.3 reduced time units were employed in the 
calculations that follow. 

The mean square displacement (MSD) is shown in Fig- 
ure[21for one-dimensional translation (the x, y, and z di- 
rections are equivalent and used as independent samples) 



TABLE I: Diffusion coefficients at each investigated temper- 
ature and density for translation and rotation as measured 
via the Einstein model and the second Legendre polynomial. 



T 




D t 


D% 


D y r 


D z 


D r (P 2 ) 


[K] 


[g/cm 3 ] 


[10~ 5 cm 2 /s] 


[us" 1 ] 


[us" 1 ] 


Ins" 1 ] 


Ins" 1 ] 


346 


1.027 


0.23 


7.6 


3.7 


3.6 


4.0 


305 


1.055 


0.057 


2.3 


1.3 


1.3 


1.0 


291 


1.065 


0.023 


0.95 


0.51 


0.58 


0.38 


275 


1.076 


0.0072 


0.42 


0.29 


0.35 


0.12 


266 


1.079 


0.0018 


0.18 


0.16 


0.24 


0.024 


260 


1.082 


0.0013 


0.14 


0.15 


0.18 


0.014 



and rotation in each of three molecular frame directions 
(see Figure . As in previous studies 0, 0> we ob- 
serve three distinct regimes for both translational and ro- 
tational displacement. Molecular motion is initially bal- 
listic with MSD proportional to (At) 2 . Following this 
initial period is a plateau corresponding to the entrap- 
ment of molecules in the cage formed by their neighbors. 
Towards the end of the plateau, the non-Gaussian pa- 
rameter attains its maximum value and thereafter the 
long-time diffusive regime begins, as evidenced by the 
emerging proportionality between the MSD and At. The 
diffusion coefficients for translation and rotation are de- 
termined from the slope of the MSD in this region and 
are listed in Table [I] An interesting feature that emerges 
from Table [I] is the relative value of D r between the 
different rotational directions. At warm temperatures, 
the rotational diffusion coefficients are ordered in ac- 
cordance with their associated moments of inertia (i.e. 
D x > D y > D z , consistent with I x < I y < I z ). This 
trend gradually reverses itself such that at 260 K, the 
diffusion coefficients obey D x < D y < D z . 

The behavior of the non-Gaussian parameter over the 
range of temperatures investigated here is shown in Fig- 
ure|5]for translation and each of the rotational directions. 
As the model OTP is cooled, At* increases and corre- 
sponds approximately to the transition from the cage to 
the long-time diffusive regime at each temperature. For 
T > 291 K, At* is approximately coincident for trans- 
lation and rotation. However, for T < 291 K, At* (T) 
increases rapidly with decreasing temperature for trans- 
lation, such that at T — 260 K, At* is an order of mag- 
nitude larger for translation than for any rotational di- 
rection. Two interesting results are shown in Figures |5] 
and [3J The first is the difference in maximum val- 
ues between the tp x and the (p y and (p z directions. The 
non-Gaussian parameter in the ip x direction reaches a 
maximum value of approximately 1.7 while a 2 in the tp v 
and (f z directions does not exceed 0.5. The second fea- 
ture of note is the non-monotonic behavior of At* (T) and 
(□^(At*) in the tp y and (p z directions. After an initial in- 
crease upon cooling down to 275 K for ip y and 291 K for 
ip z , At*(T) and a^(At*) begin to decrease as the tem- 
perature decreases. This trend continues down to the 
lowest temperature studied. In contrast, a recent study 
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FIG. 2: (color online). Mean square displacement at T = 260 (blue), 266 (green), 275 (black), 291 (orange), 305 (red), and 
346 K (magenta) for (a) translation in one dimension, (b) rotation in the tp x , (c) ip y , and (d) ip z directions. Temperature 
increases from bottom to top, and the open circle on each curve marks the time of maximum non-Gaussian behavior, At* . 



of SPC/E water showed all rotational directions to 
be qualitatively similar, with At*(T) and a^Ai*) mono- 
tonically increasing as temperature decreases. We have 
verified that these results are not an artifact of cooling 
at constant pressure as opposed to constant density, by 
performing a subset of isochoric runs and comparing the 
behavior of the non-Gaussian parameter. 

A complementary approach to the study of non- 
Gaussian behavior is based on a recently proposed al- 
ternative view of self-diffusion [2!j. Integration of the 
velocity autocorrelation function leads to the following 
expressions. 

D t = lim Uv(0)- Ar(t)) (16) 

t^oo j 

D? = lim (uj a (0)Aip a (t)) (17) 

t — >oo 



The physical interpretation of these equations is that the 
diffusion constant is a measure of the extent to which ini- 
tial velocity biases long-time displacement. Eauationsll6l 
and 1171 can be written more formally as an integral of a 
joint probability distribution of initial velocity and final 
displacement 

D t = lim [fv%AxP{v%,Ax)dv5dAx (18) 



D?= lim // to* Acp a P(u)%,A<p a )duigdA<p a (19) 

At^ooJ J 

Here P(vq, Ax) and P{oJq, Aip a ) are the joint probabil- 
ity distributions of initial velocity and eventual displace- 
ment, and the isotropy of translational motion allows the 
replacement of P(vq, Ar) with a factor of three times the 
distribution of one-dimensional displacement, P(Vq, Ax). 
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FIG. 3: (color online). The non-Gaussian parameter, 02, at T = 260 (blue), 266 (green), 275 (black), 291 (orange), 305 (red), 
and 346 K (magenta) for (a) translation in one dimension, (b) rotation in the tp x , (c) <p v , and (d) ip z directions. The time at 
which 02 attains its maximum value corresponds to At*. 



The necessary probability distributions are easily calcu- 
lated from a molecular dynamics simulation by maintain- 
ing a two-dimensional histogram of initial velocity and 
displacement after a specified time, At. 

Figure 21 shows contour plots of the joint probability 
distributions at short, intermediate, and long times at 
275 K. At short times, initial velocity and displacement 
are highly correlated, resulting in a distribution that is 
skewed along Ax — v^At and Aip a = Wq At. Interac- 
tions with other molecules soon weaken this correlation, 
and the distribution accordingly appears axisymmctric, 
with circular contours clearly visible. However, when the 
joint probability distribution is evaluated at At = At*, 
it develops a distinct "diamond distortion" . This behav- 
ior was first reported for translational motion by Shell et 
al. for a binary mixture of Lennard- Jones particles |33| , 
but has not previously been reported for rotational mo- 



tion. This distinctive shape may be reproduced by the 
combination of a Maxwell-Boltzmann distribution of ve- 
locities and a two-Gaussian distribution of particle dis- 
placements in the form of equation 1151 Random diffu- 
sive motion eventually returns the distributions to a sin- 
gle Gaussian shape at long times. We have calculated 
P(vq,Ax) and P(u)Q,Atp a ) at all investigated temper- 
atures and for several time scales in addition to those 
shown in Figure Over this range of conditions we find 
that the time of maximum "diamond distortion" corre- 
sponds approximately to At*. In addition the relative 
strength of the distortion is consistent with the values of 
a\ and a^. Namely, the contour plots become increas- 
ingly "diamond" shaped at At* as the temperature is re- 
duced with the exception of the few instances in the tp y 
and ip z directions when a decrease in temperature leads 
to decrease in the maximum value of a^. 
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FIG. 4: (color online). Contour plots of the joint probability 
distributions of initial velocity and displacement after a time 
At for (top) translation and (bottom) rotation in ip x at 275 K. 
The various values of At are marked on each plot. The <p v 
and <y3 z directions behave similarly to the ip x plots shown. 



The "diamond distortion" of the joint probability his- 
tograms evaluated at At* reveals interesting information 
about molecular motion. The circular shape of the con- 
tours at short and long times is a result of the Gaussian 
nature of the particle displacement probability distribu- 
tion, since the velocities obey the Gaussian Maxwell- 
Boltzmann distribution at all times. At intermediate 
times, characterized by At*, the tails of the measured 
probability distribution of particle displacements are sig- 
nificantly increased relative to a single Gaussian distribu- 
tion with zero mean and estimated standard deviation, 
thus causing the "diamond distortion" |(J. By comput- 
ing oi l 2 and «2 f° r various initial velocities, we found that 
non-Gaussian behavior is uniform across all initial veloc- 
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FIG. 5: (color online). Pair correlation function for the cen- 
ters of mass of mobile molecules as defined in Section IlIII at 
260 K. T-T designates the correlation between pairs of trans- 
lationally mobile molecules while R-R is the rotational coun- 
terpart with pairs of molecules simultaneously mobile in ip x , 
ip y , and tp z . T-R represents the correlation between molecules 
that are mobile concurrently in translation and all rotational 
directions. 



ities (both translational and angular) once the regime of 
ballistic particle motion has ended, after approximately 
0.4 reduced time units. This is consistent with the be- 
havior reported in for an atomic system, and confirms 
that the distinctive "diamond distortion" is entirely due 
to non-Gaussian behavior of particle displacements, and 
is unrelated to the initial velocity distribution. 

To detect the presence of spatially heterogeneous dy- 
namics we have computed the pair correlation function 
for the centers of mass of mobile molecules at 260 K, 
as shown in Figure El To compute these correlations, 
molecules are classified into mobile and immobile groups 
for translation and each rotational direction based on the 
criteria presented in Section IIIII However, the differ- 
ent values of At* for translation and the three rotational 
directions complicates this process. To allow the cal- 
culation of correlations that involve more than a single 
direction, such as correlations between translation and 
rotation or between more than a single rotational di- 
rection, molecules are classified as mobile or immobile 
based on their future displacement. This procedure be- 
gins by selecting a particle configuration and then exam- 
ining translational and angular displacements after the 
appropriate time, At*, for each direction has elapsed. 
Molecules are then classified as mobile or immobile in 
each direction, and pair correlation functions are gener- 
ated from the initial particle configuration by examining 
center of mass pairwise correlations as a function of dis- 
tance between molecules that are appropriately classified. 
This process is then repeated over many time origins to 
improve statistics. Figure shows enhanced correlations 
between molecules that are translationally mobile and 
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FIG. 6: (color online). Time dependence of the second Leg- 
endre polynomial at T — 260 (blue), 266 (green), 275 (black), 
291 (orange), 305 (red), and 346 K (magenta). The unit vec- 
tor, u, used in this calculation is shown in Figure Q Temper- 
ature decreases moving from left to right. 









between molecules that are rotationally mobile. How- 
ever, the rotational correlation is only significantly in- 
creased when it is restricted to molecules that are mobile 
in all rotational directions (i.e. molecules that are simul- 
taneously mobile in ip x , ip y , and (p z ). These enhanced 
correlations indicate that the dynamics of the Lewis and 
Wahnstrom model for OTP are spatially heterogeneous 
in both translation and rotation. In addition, the cross 
correlation between translationally and rotationally mo- 
bile molecules (i.e. molecules that are simultaneously 
mobile in translation and all rotational directions) shows 
the strong tendency for these molecules to be in close 
proximity. These results are similar to a recent compu- 
tational study of water, which found a strong simil arity 
between translational and rotational heterogeneities [15| . 

The Debye model leads to experimentally accessible 
measures of rotation, and many studies of supercooled 
liquids, including simulations, accordingly rely on it. 
Most of these examinations obtain rotational diffusion 
coefficients from the decay of the second Legendre poly- 
nomial and equation^ Figure shows the evolution of 
the ensemble-averaged second Legendre polynomial cal- 
culated for the Lewis and Wahnstrom model of OTP at 
each of the temperatures investigated. At lower temper- 
atures a two-step relaxation process is evident, as seen 
from the initial short-time decrease, which is followed by 
a long-time tail. We examined P\ through P5, but only 
show P2 because of its experimental significance. In the- 
ory, one may use any of the Legendre polynomials to 
compute D r using equation \5\ and D r = [1(1 + 1)t;] _1 , 
and should obtain the same result. However, we find 
that D r decreases as the order of the Legendre polyno- 
mial is increased (i.e D r (P\) > D r (P2) > D r (P-i) and 
so on). In addition, the Debye model predicts that the 



FIG. 7: Single-molecule rotational trajectory of the unit vec- 
tor, u (see Figure^, employed in the Debye model of rotation 
at each investigated temperature. The total time for each tra- 
jectory is the time necessary for (cosip(i)) to decay to a small 
value. 



ratio of rotational correlation times measured from the 
first and second Legendre polynomial, T1/T2, should be 
equal to 3. We find that this ratio decreases from 2.45 at 
346 K down to 1.60 at 260 K. Deviation from this the- 
oretical value in supercooled liquids is associated with 
long angular jumps |171 134| . Such behavior was shown 
to be prominent in this model at 266 K by Lewis and 
Wahnstrom [HHl. 

In addition to the Legendre polynomials, an informa- 
tive representation of the validity of the Debye model is 
the single-molecule trajectory of the vector, u, on a unit 
sphere over the time required for (cosip(t)) to decay to 
a small value [2^]. Figure {7\ shows representative trajec- 
tories for all investigated temperatures. At 346 K, u ex- 
plores uniformly the entire surface of the unit sphere in a 
manner that is consistent with the Debye approximation 
(i.e. small random steps) . Upon supercooling, the trajec- 
tory of u no longer covers the whole surface and becomes 
trapped in small regions of the sphere surface over ex- 
tended periods of time. This behavior is characteristic of 
a change in the mechanism of reorientation from consis- 
tently small random steps to well-separated and sudden 
changes of orientation. In this new mechanism molecules 
undergo librational movement for a significant amount of 
time before jumping to a new orientation. This type of 
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tions that measure rotation using the Debye model dis- 
play a significant increase in D t /D r as temperature de- 
creases. This is in accord with the traditional concept of 
translation-rotation decoupling, indicating an increase in 
the effective translational diffusion coefficient relative to 
its rotational counterpart (i.e. the system behaves effec- 
tively as if, upon cooling, molecules translate further for 
every rotation they execute). However, when D r is calcu- 
lated from the Einstein formulation, equation [5] D t /D r 
decreases as the temperature decreases, indicating an ef- 
fective increase in the rate of rotational diffusion relative 
to translation. This suggests the need for a critical re- 
examination of our current understanding of translation- 
rotation decoupling in supercooled liquids, especially in 
light of its dependence on the Debye model. 



FIG. 8: Temperature dependence of the logarithm of the ratio 
of translational to rotational diffusion coefficients normalized 
by the corresponding value at high temperatures. Experimen- 
tal data are from fla.ll4|. 



movement is clearly visible in the trajectory of u at 266 
and 260 K. Figure vividly shows a breakdown of Debye 
behavior. 

One of the distinctive aspects of supercooled liquid be- 
havior is the decoupling of translational diffusion from 
viscosity (breakdown of the Stokes-Einstcin equation) 
and of rotational diffusion from translational diffusion . 
It has been the subject of many experimental (for OTP 
see 0,0>HSEi3) and theoretical studies [Tlllr2ll3^ . l38| . 
Above the melting temperature D t and D r are propor- 
tional to T/rj. This is in accordance with the Stokes- 
Einstein (SE) and Debye-Stokes-Einstein (DSE) equa- 
tions (see equations ^ and |5] respectively) . Experiments 
on OTP and other deeply supercooled fragile liquids show 
that the SE relation seriously underpredicts D t whereas 
the DSE relation remains substantially valid down to 

A convenient metric of translation-rotation decoupling 
is the quotient D t /D r normalized by the same quantity 
at high temperature |l2l (one may equivalently use 
D t Ti in place of D t /D r ). Figure |H1 shows this quantity as 
a function of temperature, using diffusion coefficients ob- 
tained at 346 K as the high-temperature reference values. 
Included in this figure are rotational diffusion coefficients 
based on the Debye model obtained from (P2[cosip(t)]) ) 
as well as the three rotational coefficients resulting from 
the Einstein formulation. Figure |S| also includes ex- 
perimental rotational correlation times calculated from 
deuteron spin alignment ( H-NMR) experiments using 
equation 01 with I = 2 |l3ll39l|. and experimental transla- 
tional diffusion coefficients from a recent study by Mapes 
et al. 0]. 

Figure |S] reveals a significant result. Upon entering 
the deeply supercooled regime, experiments and simula- 



V. CONCLUDING REMARKS 

The diversity and complexity of dynamic phenomena 
present in supercooled liquids are a major challenge to 
a comprehensive understanding of this important class 
of condensed-phase systems. Computer simulation is an 
important tool that provides insight into the details of 
molecular motion in a manner that is not currently pos- 
sible in experiments. To further the understanding of 
diffusion in supercooled liquids, we have studied a model 
of the canonical fragile glass-former ort/10-terphcnyl. In 
particular, we have extended the formalism and tech- 
niques developed for studying dynamic heterogeneity in 
translational motion [6J, |29J to a molecular system with 
rotational degrees of freedom. These methods revealed 
spatially heterogeneous dynamics in translation and rota- 
tion with a strong spatial correlation between the trans- 
lationally and rotationally mobile molecules. The com- 
monly used Debye model of rotation was shown to break 
down at deeply supercooled temperatures, as the mecha- 
nism for molecular reorientation begins to incorporate 
large angular jumps. When the Einstein formulation 
of rotational motion was used to examine translation- 
rotation decoupling, the analysis showed a trend oppo- 
site to that observed when using the Debye model to 
quantify rotational diffusion. Specifically, the effective 
rate of rotational motion appears to be enhanced rela- 
tive to translation. This result, coupled with the concur- 
rent breakdown of the Debye model, calls into question 
conventional interpretations of the relationship between 
translational and rotational motion in deeply supercooled 
liquids. 

Models that explain translation-rotation decoupling 
are based on t he p icture provided by dynamic hetero- 
geneity [H [Hli3 and rely on the Debye model to de- 
scribe rotation. By assuming the presence of regions of 
fast and slow dynamics, translation-rotation decoupling 
emerges as temperature falls below some critical value 
(e.g. T < 1.2T ff ) as a consequence of the different ways 
in which translational and rotational motion are averaged 
in regions of slow and fast dynamics [Tlj . In this view, 
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the SE and DSE equations are assumed to be obeyed 
locally in both the slow and fast regions of the dynami- 
cally heterogeneous liquid 0. It can then be shown that 
the effective translational diffusion coefficient is given ap- 
proximately by (D%+D{)/2, where the superscripts s and 
/ denote the slow and fast regions, respectively, and it 
is assumed in this approximate calculation that the slow 
and fast regions each account for half of the system's vol- 
ume 0. Using similar arguments, the rotational correla- 
tion time is given by (r/ 5 + t, )/2. Thus, the translational 
diffusion coefficient is determined by the dynamics of the 
fast regions, whereas the rotational correlation time is 
determined by the dynamics of the slow regions [lj. A 
corresponding microscopic interpretation of the results 
shown in Figure [S] when the Einstein formalism of rota- 
tion is used has yet to be developed. In particular, an 
understanding of how heterogeneity affects averaging in 
such a way as to produce an effective enhancement of 
rotation upon cooling needs to be developed. 

An interesting question arising from our work is the 
origin of the non-monotonic behavior of At* and a^iAt*) 
in the ip v and (/^-directions [see Figures EJc) and|3Jd)]. 
It is possible that this behavior may be caused by the on- 
set of orientational hopping in the ip v and ^-directions. 
Future work will focus on the onset of orientational hop- 
ping in the various directions, and in particular will test 
the eventual appearance of hopping in the (^-direction. 

This and previous studies of the Lewis and Wahnstrom 
model for OTP suggest that its behavior differs from that 
of real OTP at supercooled temperatures. Its primary 
fault is the inaccurate prediction of diffusion coefficients. 
Our study and others |2J, [Mj report diffusion coefficients 
that are three orders of magnitude larger for translation 
and seven orders of magnitude larger for rotation than 
experiments near 260 K indicate. This may in part re- 
sult from fitting the Lennard- Jones interaction parame- 
ters to experimental values for the translational diffusion 
coefficient and molar volume at 400 K [24| as opposed to 
a lower temperature. This change could be made with 
relative ease, but raises the question of which molecular 
features of OTP contribute most to its glass-forming abil- 
ity. OTP is known to interact with short-range van der 
Waals forces [4JJ, and the molecular structure exhibits 
some internal torsioning. These features have been in- 
corporated into other models for OTP, including an 18 



Lennard-Jones site, non-rigid molecule |42( and a model 
that uses fully atomistic force field methods to describe 
the interactions 43] . While more accurate, particularly 
in the fully atomistic case, these models make it compu- 
tationally challenging to use the system sizes and simula- 
tion times necessary at deeply supercooled temperatures. 
The prominence of OTP as one of the most extensively 
studied fragile glass-formers attests to the importance of 
developing an accurate model that captures the salient 
features of real OTP but is simple enough for use in sim- 
ulation studies at supercooled temperatures. The Lewis 
and Wahnstrom model is a first step towards this goal, 
but improvements are warranted. 

Our present analysis suggests further questions regard- 
ing the nature of dynamic heterogeneity. An important 
open question is how regions of high mobility emerge in 
the liquid. It is not yet understood what local proper- 
ties of the liquid cause the molecules in these domains 
to have a high mobility. An interesting method to ex- 
plore this question was introduced by Widmer-Cooper et 
al. [44l l45l |46| . This technique involves running sepa- 
rate simulations from a single starting configuration. At 
the beginning of each run the momenta of each parti- 
cle are randomly chosen from the appropriate Maxwell- 
Boltzmann distribution. Their results indicate that a 
particle's local environment, and not its initial velocity, 
has a strong effect on its propensity for motion. This sug- 
gests the existence of structural features that influence 
particle mobility and engender dynamic heterogeneity. 
It would be interesting to extend a study of this nature 
to include rotational degrees of freedom. The Einstein 
rotational formalism used in this paper lends itself well 
for such a study. Knowledge of the origin of dynamic 
heterogeneity would be a significant advance in the un- 
derstanding of supercooled liquids. 
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